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We show that the surface of an s-wave superconductor decorated with a two-dimensional lattice 
of magnetic impurities can exhibit chiral topological superconductivity. If impurities order ferro- 
magnetically and the superconducting surface supports a sufficiently strong Rashba-type spin-orbit 
coupling, Shiba sub-gap states at impurity locations can hybridize into Bogoliubov bands with non¬ 
vanishing, sometimes large, Chern number C. This topological superconductor supports C chiral 
Majorana edge modes. We construct phase diagrams for model two-dimensional superconductors, 
accessing the dilute and dense magnetic impurity limits analytically and the intermediate regime 
numerically. To address potential experimental systems, we identify stable configurations of fer¬ 
romagnetic iron atoms on the Pb (111) surface and conclude that ferromagnetic adatoms on Pb 
surfaces can provide a versatile platform for two-dimensional topological superconductivity. 


The chiral p-wave superconductor in two dimensions 
(2D) and the closely related fractional quantum Hall 
Pfaffian state at v = 5/2 are the archetypal examples 
of topologically ordered states of matter that support 
non-Abelian anyonic excitations. [1-3] The theoretical 
exploration of these states has shaped our understanding 
of topological order and is foundational for the distant 
goal of building a topological quantum computer. [4- 
6] In contrast to the v = 5/2 fractional quantum Hall 
state, however, to date chiral p-wave superconductiv¬ 
ity has not been confirmed in any experimental system. 
The most prominent candidate system, superconducting 
Sr 2 Ru 04 , [7] has been the subject of an ongoing debate 
about its actual paring state. [8-11] 

Chiral superconductors break time-reversal symmetry 
(TRS). This hinders the formation of Cooper pairs, since 
orbital (and possibly paramagnetic) pair-breaking effects 
can come into play. Depairing is also the main hurdle for 
realizing a line of proposals in which layered heterostruc¬ 
tures involving ferromagnets and s -wave superconductors 
are used to build an artificial 2D p-wave superconduc¬ 
tor. [12-14] The guiding principle for these proposals is 
to design a band structure with a single normal-state 
Fermi surface with no spin-degeneracy. If Rashba spin- 
orbit coupling is present so that the states on this Fermi 
surface are not fully spin-polarized, and if s-wave super¬ 
conductivity is proximity-induced in such a system, the 
effective pairing near the single Fermi surface is equiva¬ 
lent to that of a chiral p-wave superconductor. 

In one dimension (ID), based on the principle of com¬ 
bining spin-orbit coupling and externally applied mag¬ 
netic fields, various groups have proposed engineering 
artificial realizations of p-wave superconductors. [15-21]. 
An experimental realization of these proposals employ¬ 
ing semiconductor nanowires with strong spin-orbit cou¬ 
pling [22] has reported Majorana fermion signatures. In 
this setup, the externally applied magnetic field has to 
be rather small (to avoid suppressing superconductivity). 



FIG. 1. Concept of building a topological p-wave supercon¬ 
ductor from magnetic impurity bound states on the surface 
of an s-wave superconductor like Pb with strong spin-orbit 
coupling. A very thin layer of magnetic adatoms, such as Fe 
or Co, is deposited on the surface of Pb, and forms islands 
around which chiral Majorana edge modes are present (yel¬ 
low). Depending on the density of the magnetic atoms, one or 
many Majorana modes may appear. A one-dimensional chain 
of magnetic atoms may support localized Majorana modes at 
its end. These modes can be detected by tunneling tech¬ 
niques. 


As the phase-space for the existence of a topological su¬ 
perconductor is controlled by the Zeeman gap [23] , these 
systems require a delicate balance of the parameters in¬ 
volved (spin-orbit coupling, magnetic field and chemical 
potential) in order to create the topological superconduc¬ 
tor. 

Recently, a ID topological superconductor was real¬ 
ized in a system that is quite distinct but employs sim¬ 
ilar microscopic ingredients - spin-orbit coupling, ferro¬ 
magnetism, and s-wave superconductivity. [24] A chain 
of magnetic Fe atoms is deposited on the surface of an s- 
wave superconductor with strong spin-orbit interactions. 
The Fe chain is ferromagnetically ordered [24] with a 
large magnetic moment, which takes the role of the mag¬ 
netic field in the nanowire experiments. Unlike previous 
proposals, this “magnetic field” is mostly localized on the 
Fe chain, with small leakage outside. Superconductivity 
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is not destroyed along the chain. In this setup, the energy 
scale of the exchange coupling of the Fe atoms is typically 
much larger than that of the Rashba spin-orbit coupling 
and the superconducting pairing. The ferromagnetically 
ordered Fe atoms induce localized Shiba states within the 
gap of the superconductor. [25-27] The hybridization of 
these states forms the band structure of a ID p-wave 
superconductor that supports Majorana end states. [28] 
Because the Fe bands are fully spin split, no additional 
control over the chemical potential is necessary. A sim¬ 
ilar scenario applies when the Fe orbitals are magnetic 
but itinerant. [23] 

In this Letter we point out that this strategy can also 
be successful in 2D. Magnetic adatoms on the surface of 
a superconductor with strong spin-orbit coupling, when 
arranged in a 2D lattice, [29] can yield a 2D topologi¬ 
cal chiral p-wave superconductor whose chiral Majorana 
edge modes can be observed in Scanning Tunneling Mi¬ 
croscope (STM) measurements (see Fig. 1). To shed 
light on the rich range of possibilities, we analyze the 
topological properties of a system with dense local mo¬ 
ments that are exchange coupled to a model 2D supercon¬ 
ductor, demonstrating that topological superconductors 
with higher Chern numbers, and consequently multiple 
chiral Majorana edge channels, can easily occur. We are 
also able to analyze the model’s dilute magnetic impu¬ 
rity limit analytically and obtain numerical topological 
phase diagrams for intermediate impurity concentrations. 
Based on density-functional-theory (DFT) calculations, 
we further propose realizing 2D topological chiral p-wave 
superconductors experimentally by depositing transition 
metal adatoms on superconducting Pb. The type of mag¬ 
netic ion can be varied to access different strengths of 
the magnetic moment. In the case of Fe adatoms on a 
Pb (111) surface, we show that strong magnetic order in 
general leads to an odd number of 2D Fermi surface seg¬ 
ments. As a consequence the proximity-induced super¬ 
conducting phases can have nonzero Chern numbers and 
chiral Majorana edge modes. Below we first present our 
model system results which bring out a number of generic 
features of superconductor surfaces with ferromagneti¬ 
cally ordered magnetic adatoms, and then we summarize 
our DFT results by focusing on the specific case of Fe 
adatoms on a Pb (111) surface. 

To render the analytical calculations tractable, we 
consider a Hamiltonian that models only the surface 
layer of a bulk 3D 5-wave superconductor on which 
the 5-wave superconducting order parameter A is in¬ 
duced from the bulk. On this superconducting layer, we 
model the magnetic impurities as classical spins whose 
only interaction with the electrons in the superconduc¬ 
tor is through Zeeman-like couplings. [25-27] Employing 
a tight-binding description on a 2D square lattice A that 
is spanned by the primitive lattice vectors e\ and & 2 , we 


consider the mean-field Hamiltonian 

H = ^ ^ 1 1 (cJ,Cr+ei "b C r C r+e 2 ) — ~^ C r C r T A 4,t C r*,t 
re A 

+ia (cl(T2Cr+& t - 4<TiC r+ e 2 ) + h.C.] 

+ j E ci<7 3 C r . 

re A* 

(i) 

Here, cj. = (4 4 |) * s a s Pi nor °f the creation operators 

for electrons at site r with spin and cq, i = 1,2,3, 
are the three Pauli matrices. We denote by t the near¬ 
est neighbor hopping integral in the superconductor, /i 
the chemical potential, and a the strength of the Rashba 
spin-orbit coupling. Classical magnetic moments ferro¬ 
magnetically aligned normal to the plane in the 03 di¬ 
rection are positioned on a sublattice A* of A and are 
exchange-coupled to electrons via the term proportional 
to J. Drawing experience from the ID situation and the 
ab-initio calculations presented below, the physically rel¬ 
evant hierarchy of energy scales that we consider here is 
given by 

£ » J » a ~ A. (2) 

Dense impurity limit — As a warmup, it is instructive 
to consider the simplest situation in which each lattice 
site is coupled to a magnetic moment, i.e., A = A*. This 
limit is representative of system with self-assembled is¬ 
lands of magnetic adatoms. Its consideration allows us to 
highlight the difference between the regime Eq. (2), and 
that of small J that has been studied previously. [12-14] 
In particular, it is possible to access a phase with Chern 
number 2 in the large J limit. For the case A = A* and 
a = 0, the Hamiltonian (1) has gapless lines in momen¬ 
tum space defined by 

e k = ±yj 2 - a 2 , (3) 

where = 2£(cos Aq+cos A^)—/i. Adding spin-orbit cou¬ 
pling a / 0 will generically lift these degeneracy lines and 
gap the spectrum except at the four inversion-symmetric 
momenta [kij = (1 — i, 1 — j) 7r/2 with i, j = d=] at which 
the spin-orbit coupling vanishes and at which phase tran¬ 
sitions between different Chern number phases occur. At 
each of these momenta, the condition (3) is met for two 
values of the chemical potential 

a = —2 (i + j)t + AvO 2 - A 2 , i,j, A = ±. (4) 

Around each of these gapless points in the three- 
dimensional fi-k parameter space, we can reduce the 
Hamiltonian to an effective two-band model and expand 
it to linear order in the deviations 5k from and the 
deviations 5fi from yielding 


= bk^cri - i5ki<j 2 ) + Ay 1 - 5fur 3 . 

( 5 ) 
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FIG. 2. Phase diagrams for Hamiltonian (1), in terms of 
Chern numbers [(a) and (b)] and superconducting gaps [(c) 
and (d)] as a function of the chemical potential /x and the 
strength of the magnetic moments J. The panels (a) and (c) 
[(b) and (d)] correspond to the case of one magnetic adatom 
every 2x2 (3x3) lattice sites. In these plots, A =? 0.06t and 
a = 0.11; /x has been shifted such that y — 0 implies the 
chemical potential lying in the center between the two spin- 
split band bottoms. 


potentials, more than two Fermi surfaces can exist in the 
reduced Brillouin zone defined by A* as the Fermi wave¬ 
lengths are significantly smaller than the lattice spacing 
of A*, as a consequence higher Chern numbers can occur 
[e.g. 8 in Fig. 2(a) and 15 in Fig. 2(b)] but the trade¬ 
off is an overall smaller induced gap; iii) the phases with 
different Chern numbers are generally separated by lines 
(in the 2D parameter space) defined by conditions simi¬ 
lar to Eq. (4), and across each specific line the change of 
Chern numbers is a constant determined essentially the 
same way as in our preceding analysis. 

Dilute impurity limit — To better understand the di¬ 
lute impurity limit, we complement our results on the 
lattice by a calculation in which we treat the underlying 
2D superconductor in the continuum limit and consider 
sparsely distributed Shiba impurities that are arranged 
in a square lattice on top of it. This allows us to derive an 
effective two-band model for the hybridizing Shiba states. 
This effective Hamiltonian represents a chiral p-wave su¬ 
perconductor in the appropriate parameter regime. 

The strategy of our derivation is inspired by the cal¬ 
culation for ID Shiba chains of Ref. [28] (for details, 
see [30]). We start from a 4 x 4 Bogoliubov-de-Gennes 
(BdG) Hamiltonian 


Here, the Pauli matrices act on a subspace defined by the 
two bands that satisfy condition (4). In the y-k space, 
Hamiltonian (5) represents a source of Berry flux that 
corresponds to a unit topological charge i x j x A = ±1. 
As the chemical potential ramps through the gapless 
point at the Chern number changes by i x j x A. Let 

us assume 2 1 > VJ 2 — A 2 , so that the critical chemical 
potentials are in the order /x_ 5 _ ? _ < /x_ 5 _ 5+ < /x_ 5+5 _ = 

M+< M-,+,+ = M+< D +yf,- < M+,+Hr T^- 
ially, for y < /i-,-,- the Chern number is zero. By 
increasing /x, the system thus exhibits phases with the 
Chern numbers 

C= -1, 0, +2, 0, -1, (6) 

where the phase transitions occur at the respective /i^y. 
Hence, with homogeneous magnetization, the supercon¬ 
ductor may already exhibit Chern number equal to 2 . 
In cases where magnetic impurities are spaced more 
sparsely, i.e., if A* is a sublattice of A, even higher Chern 
numbers can be obtained. 

We solved Hamiltonian ( 1 ) numerically for the case of 
one magnetic impurity every 2 x 2 and 3x3 plaquettes of 
the square lattice, and show the phase diagrams in Fig. 2. 
From the phase diagrams, we can read off three general 
features: i) at small chemical potentials around the band 
bottom, where the Fermi wavelengths are larger than or 
comparable to the lattice spacing of A*, the sequence of 
Chern numbers (—1,0, + 2 ) always occurs when J ft is not 
too large - this universal feature corresponds to the dense 
limit that we have discussed above; ii) at larger chemical 


=(a- 4) (7a) 

that acts on 4-spinor valued wave functions 4/(r) = 

W’f>V’l,V4—V’j) T ( r )> r e r2 > with = k 2 /(2m) - fi + 
<y(ki<j 2 — k 2 (Ti). In addition, the magnetic impurities are 
represented by the Hamiltonian 

H imp = -J Yi 6(r-r*)S 3 , (7b) 

r*eA* 

where S 3 = 03 ® To, with T 3 the identity matrix acting on 
particle-hole space. If we restrict the wave-function T(r) 
to the locations r* G A* of the impurities, we obtain the 
self-consistency equation 

(E - H* ag ) _1 (- JS 3 )V (q) = $(g) (8) 

for the Fourier transforms ^(q) = J^ r * eA * e ~ iqr T(r*), 
where the momentum q G [0, 2tt) 2 now belongs to the A* 
Brillouin zone and we have set the impurity spacing to 
unity. 

We need two more steps to reduce Eq. ( 8 ) to an ef¬ 
fective two-band model for the Shiba states, assuming 
they are deep in the superconducting gap and dilute 
compared with the Fermi wavelength. First, the left- 
hand side is expanded to linear order in the energy E 
to cast the equation in the form of the time-independent 
Schroedinger equation. Second, we project the effective 
Hamiltonian into the eigenstates of an isolated Shiba im¬ 
purity on every site, given by 4/ + = (1,0, l,0) T /y/2 and 
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FIG. 3. Parity of the Chern number for the BdG band 
structure of an s-wave superconductor decorated with dilute 
Shiba impurities, following Eq. (10). Here, a a* is the lattice 
spacing of the Shiba impurities, Af is the Fermi wavelength 
in the limit a 0 and the Hamiltonian (9) was evaluated 
with £ — 30 Af- 


= (0,1,0,-1) T /V2 | 

band Hamiltonian 


We obtain the effective two- 


q V + d 0:Q 


[(1 -Tj- d 0l q) cr 3 - d 2 , q cr 1 + di jQ <7 2 ] 


where 


d ^ = i E e_i9 ' r */ + (kl), 


(9a) 

(9b) 


r *-^0 


di, q =^E e ~ iq ' r '&f-(\ r *\)’ ( 9c ) 


r*^ 0 


are defined in terms of the functions 


f±(r) 


a-r/€ 


Vp 


^ 61 


V(±1) C1 


-A)/2. 


'k'f 1 cos ( kp 


(4 A) r - J) , 


(9d) 

In the equation above, £ is the coherence length of 
the 2D superconductor (without the magnetic impuri¬ 
ties), vp = + 2/i/m is its Fermi velocity, = 

^/2m/i + <a 2 m 2 =p rna are the Fermi wave-vectors for its 
two spin-split bands, and r] = J[k^ + k^]/Avp is a 
dimensionless parameter. In the “deep Shiba limit” in 
which the projection in the states and 4/_ is justi¬ 
fied, we have 77 rsj | # 

The Hamiltonian (9) represents the effective super¬ 
conductor formed by the Shiba bound states within the 
gap of the underlying s-wave superconductor. Similar to 
the case of the effective two-band model in Eq. (5), this 
Hamiltonian can have nodal points in k-fi -space at which 
the Chern number changes. However, unlike in Eq. (5), 
the nodes can occur at any point in the Brillouin zone, 
making an analytic treatment intractable. In addition, 


the validity of Hamiltonian (9) requires a self-consistency 
that permits the low-energy expansion of Eq. ( 8 ). There¬ 
fore, instead of computing the Chern numbers in an ex¬ 
tended parameter space, we focus on information that 
can be obtained at special points of the Brillouin zone 
at infinitesimal energy. To that end, observe that Hamil¬ 
tonian (9) has C 4 rotational symmetry. Thus, any gap 
closing at points other than the (^-symmetric momenta 
k = (0, 0) and k = ( 7 r, 7 r) changes the Chern number by 
an even integer due to the symmetry-imposed multiplic¬ 
ity of the nodal points. By expanding the Hamiltonian 
into Dirac form around the C 4 -symmetric momenta, we 
obtain the expression [30] 

(- 1 ) C = sgn [d 0 , ( 0 , 0 )] sgn [d 0)(7ri7r) ] (rj = 1 ) ( 10 ) 

for the parity of the Chern number. The numerical eval¬ 
uation of this equation is shown in Fig. (3) in the form 
of a phase diagram. 

We have also performed calculations for magnetic or¬ 
ders other than simple ferromagnetism. In particular, 
the case where the magnetic configurations corresponds 
to 2D helices is related to previous studies on ID he¬ 
lices. [31-40] We obtained criteria for such a system to 
be fully gapped by proximity effect, and found that the 
fully-gapped superconducting phases can be generically 
topologically nontrivial. The results and phase diagrams 
are presented in the supplementary material. [30] 

In closing, we complement our simple model consider¬ 
ations with a specific material proposal to realize a su¬ 
perconductor with nonzero Chern numbers. For that, 
we consider transition metal atoms, in particular Fe, de¬ 
posited on the (111) surface of Pb, a strong type-II su¬ 
perconductor. The same combination of materials, but 
a different surface of Pb, was used in the experimental 
realization for the ID p-wave superconductor. [24] The 
Pb atoms on the (111) surface form a triangular lattice. 
Through ab-initio calculations that assume a ferromag¬ 
netic alignment of the Fe magnetic moments and include 
spin-orbit coupling, we compared the relaxation energies 
for various densities and arrangements of Fe adatoms on 
the Pb (111) surface, and found that a deposition with 
one Fe atom in each triangular plaquette is particularly 
favorable (for details, see [30]). In this case, the Fe atoms 
form a honeycomb lattice in which the atoms sit at differ¬ 
ent heights in each sublattice [see Fig. 4(a)]. We further 
performed DFT calculations of the electronic structure. 
The resulting (not spin-degenerate) Fermi surface and 
band structure restricted to the Fe d-orbitals is shown in 
Fig. 4(b) and (c). Critically, we find an odd number of 
Fermi surfaces. The Chern number of the corresponding 
BdG Hamiltonian is indeed nonzero over a large range of 
the chemical potential [see Fig. 4(d)]. 

In conclusion, we have proposed a versatile platform 
for realizing chiral superconductors in 2 D. We have ob¬ 
tained analytically the topological phase diagram (Chern 
number and gaps of the superconductor) of the dilute and 
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FIG. 4. Electronic structure from DFT calculations for Fe 
adatoms on Pb (111) surface, (a) The atomic configuration: 
Fe adatoms form a buckled a honeycomb lattice on the Pb 
(111) surface, (b) Fermi surface structure of the Fe d-orbitals 
for the Fermi energy indicated by the dashed lines in panel 

(c) and (d) (see [30] for cases with different Fermi energies). 
A single Fermi surface is found around the T point, (c) DFT 
band structure of the Fe d-orbitals with spin-orbit coupling. 

(d) Chern number of the BdG bands when a small s-wave 
superconductivity (A = O.OleU in this calculation) is induced 
in the system. 


dense limit, and numerically evaluated the phase diagram 
in the intermediate regime. We then showed through 
a more realistic ab-initio calculation that ferromagnet- 
ically ordered Fe atoms on the (111) surface of Pb in 
the dense limit could give rise to a chiral superconduc¬ 
tor. The presence of a 2D chiral superconductor could 
be established experimentally by tunneling into the chi¬ 
ral Majorana modes, whose number is equal to the Chern 
number of the phase, and which would take place only 
on the edge of a 2D thin island of Fe on the surface of 
Pb. 

Note During the completion of this work, a similar 
proposal for realizing chiral superconductivity appeared 
on the arxiv. [41] Where the present manuscript overlaps 
with [41], in the dilute limit, our calculation agrees with 
the one presented in [41]. Our paper also further presents 
the non-dilute limit and the ab-initio results for realistic 
materials. 
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Derivation of the effective Hamiltonian in the dilute Shiba limit 

General strategy 

We want a general formalism to find the low energy excitations and effective Hamiltonian of the following Hamil¬ 
tonian 


H = H 0 + H lt H^Y.VnSir-rn), (11) 

nGc 

where Hq is the original Hamiltonian which is gapped (« A) around zero energy, V n are matrices associated with the 
inner degrees of freedom (spin, particle-hole ect.) which can induce in-gap states, and r is implicitly d-dimensional 
(also n, m, &, g, R, n 5 below). 

We start with the Schrodinger equation 


(H 0 + H 1 )V = E9, (E < A). 


( 12 ) 


It follows that 


G 0 (E)H 1 $ = 


(13) 


where Gq{E) = (E — iTo ) -1 is the Green function. 

Because Hi is composed of delta-functions for a small set c, GqHi is nonzero only in the columns corresponding to 
c, thus Eq. 13 is equivalent to 


E G o ( E ; r , r n )V n ip(r n ) = ip(r). (14) 

nGc 

In the simplest case, if c contains only one single point, labeled by 0, then Eq. 14 implies 


G 0 (E;r 0 ,r 0 )V 0 ip(r 0 ) = ip(r 0 ). 
The energy of the excitations can be obtained by solving 


(15) 


Det[l - G 0 (E;r 0 ,r 0 )Vo] = 0 , 

Go(E-r 0 ,r 0 ) = Go(E-,0)= j dk G 0 (E;k). 


(16) 

( 17 ) 


In more complicated cases, the following equation, again implied by Eq. 14, can serve as the starting point to extract 
an effective Hamiltonian 


Vra G c : E (18) 

nGc 

G 0 (E;r m ,r n )=G 0 (E-r m -r n ) = j dk G 0 (E; k)e ik ^~ r "\ (19) 

In addition if Vn : V n = Vq and r n = nR (n G with d the dimension), Eq. (18) can be transformed to k- space: 


where 


V’O?) = E e ^"'V’O'™) (<? e [—7r/J2, tt/R)) 

m 

= J2 e ~ iqrm Go(E-,r m - r n )Voip(r n ) 

m,n 

= E/ dfc G o(^;fc)e i ( fe -' ?)(rm - r - ) Foe- i9r '*V’(r-n) 


= J dk G 0 {E- k) 

= J dk G 0 (E- k) 

= G 0 (E;q)V 0 ^(q), 


E 

_ m' 


0 i{k—q)m'R 


V 0 


J2e~ iqr ^(r n ) 


(5(fc — q — 2ribTr/R) 


Vo1>(q) 


G 0 {E-, q) := ^ G 0 (£ ; ; 9 + 2n b n/R). 


( 20 ) 

( 21 ) 

( 22 ) 

(23) 

(24) 

(25) 

(26) 


ni) can be interpreted as the band index when the /c-space is folded into the Brillouin zone defined by [—ir/R,ir/R). 


Shiba lattice in the 2D Rashba superconductor 


We start with the BdG Hamiltonian 


H BdG (k) = 



fr? 

6c = 2ffl* ^ M ~b Oi{k x (7y ky(J x ) , 

A/c — ^s^"0 H - fipfax&y kyO'x')' 


First we make the Hamiltonian dimensionless by defining 

kft = y2m*/u/fi 2 , k = k/k li , 
la = k? I (2m* a), a = l/(l a kp), 
l P = h 2 /(2m*5 p ), A p = 1 /(Ipkp), 
H = H/n, A g = A s /fi, 


and then omit the tildes, the bulk Hamiltonian becomes 

HB dG (k) = T Z ® [(k 2 - 1)<t 0 + a(k x a y - k y a x )] 

“t - T x ® [A s cto 4“ A p (k x a y &y^x)]- 


(27) 

(28) 
(29) 


(30) 

(31) 

(32) 

(33) 


(34) 


In this dimensionless Hamiltonian, energies are measured in units of fi (assumed to be large compared with other 
energy scales), lengths (wave vectors) are measured in units of the Fermi wavelength (wave vector). 
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Next we go to the basis where the normal state Hamiltonian is diagonalized 

H BdG (k ) = U{vk) ] H BdG {k)U{<p k ) 

= 7+0 [( k 2 - 1)cto + aka z \ + 7+0 [A s ct 0 + A p ka z ] 
= \E+(k)r z + A + {k)r x ] © [E_{k)r z + A_(fc)r x ], 


where 


U(<p k ) To © ^ \i e iVk ie iVk J > 

-E±(/c) = k 2 — 1 =b <xk , 

A±(fc) = A s =b A p k. 


According to the general formalism (Sec. ), we need to evaluate 

G 0 (E,r) = J dk G 0 (E,k)e ik r 


r27r dip k 


Jo ^tt Jo 




C ^ U{v' k )nE,r cos v' k )U{v' k ) ] 


R(ip, 




where 


r (e,x)= [ 
Jo 

R(<Pr) =Tq® 

To evaluate T(E,x), we first perform the integral 


00 k dk 

e ikx [E - H BdG {k)]-\ 


27T 
1 0 

0 e ^ 7 


r°° lu fiu 

T b=± (E,x) = j [ —e to [£ ~ E b {k)r z - A^t*]” 1 

~ Pb J dE b e ^ ( F b) +E b / VF]x E + E b r z + A b r x 


E 2 -E 2 -A 2 ’ 


(35) 

(36) 

(37) 


(38) 

(39) 

(40) 


(41) 

(42) 

(43) 

(44) 

(45) 

(46) 

(47) 


where vp = 2 /.:q ( ko = */l + o; 2 /4) is the Fermi velocity (the same for ± bands), /.:^ ±) = ko © a/2 is the Fermi wave 
vector for +/—band, p b = kp^/2nvp is the density of states for 6-band at the Fermi energy, and A;, = A b (kp^) is 
the pairing gap at the Fermi energy. Following Pientka et al. [42], we have 


J dE b e iEbX ' VF 


1 


E? + a 2 a 


— ^ e -a\x\/v F 


E h 


UJ 


D 


E% + a 2 E%+u,* d 


U D 


uj d — a* 


sgn(x)(e 


-a\x\/v F _ p -uj d \x\/ v F 


), 


(48) 


(49) 


where a = a/ A% — E 2 is a short-hand notation, and u> B is the Debye frequency which will be sent to +oo in the end. 
Then 


r& — 706 + JlbTx + 736 T z , 

( 7o+ + 73+ 7i+ 0 \ 

0 7o- + 73- 0 7 i_ 

7i+ 0 7o+ - 73+ 0 

0 71— 0 7o- — 73-/ 


r = r+ © r_ = 


(50) 

(51) 


7i- 
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where 


E 


7o b = ~npb 

7l6 = -npb 


ik ^ x -^/ A 2 - E 2 \ x \/ v F 
^ b pik^ x-^/A 2 -E 2 \x\/v f 


736 = -7T/>6 ~2 


U D 


UJ- 


D 


+ E*-A 2 h 


sgn(x)e ik F )x (e-^ A t- E 2 W /vF - e ~^n\ x \/v F y 


Next we find 


U(p k )TU(p k ) ] = ^ 


1 


A+ -ie~^ k A_ C+ —ie~ iVk C-\ 
ie i<fik A- A+ ie iv>k C- C+ 

C+ —ie~ i(pk C- B+ —ie~ i<fk B- 
\ie iVk C- C+ ie iipk B_ B + J 


where 


A ± = (7o+ + 73+) ± (7o- + 73-), 
B± = ( 70 + — 73+) ± ( 70 - - 73-), 
C± = 7i+ ± 7i-• 


The integral Eq. (43) breaks down to the following basic one 

/>7r/2 

h(z)m ' 


(z) = i l f dp e <('v+*co6 V ) = Q) ±1 ) 

J—ir /2 


r —7r/2 
ptt /2 

= i l / dpe"* 

J — 7T / 2 

+ OO 


+ OO 


J] 


n =—00 
■7t/2 


(Jacobi-Anger expansion) 


~r^ r 77 /- 

= 51 * ,+ ".7 n (;s) / d^e i( ' +n)v 

n =—00 — 7 r /2 


+00 


= ttJ-i(z) + 2i 51 2to+1 , 


m =—00 


where J n (z) is the the n-th order Bessel function of the first kind. Expressed in terms of Ii(z ), we have 

J o ~2^~ ellVk lob{ cos fk) = -y y A 2^_ ^2 [6(*&) + ( _1 )*6(-4)] » 

^-i l e av >'n b (co8<p k ) = -f ^j= [Ii(z b ) + (-l) 1 /,^*)] , 

J yy i z e dVfe 7 3b (cos^ fc ) 

= -y ^ + ^L A 2 IW - (-l)'*(-4) - 6(4) + (-l) , / I (-4*)] - 

2:6 = + i-^/ — Er/vp, z' b = kp\ + iLo^r/vp- 

First we look at the limit r —» 0, such that z b , z' h —> 0. In this case 

Io(z = 0) = 7t, I±i(z = 0) = ±2 i. 

We find 

Go(E,r -+ 0) = i?(</? r )(4ro ® cr 0 + C r x ® a 0 )R(p7 


(52) 

(53) 

(54) 


(55) 

(56) 

(57) 

(58) 

(59) 

(60) 
(61) 
(62) 

(63) 

(64) 

(65) 

( 66 ) 

(67) 

( 68 ) 


(69) 





















11 


= A To ® cr 0 + C T x 0 ( 7 0 , 



d(fk 

2 i r 


(7o+ +7o-) 


r—>-0 




d^k 

2 i r 


( 71 + +7i-) 


r—>■() 


it 

^ 7T ^ pb^b 

[ ~2>^^TW 


( 70 ) 

( 71 ) 

( 72 ) 


For r > 0 , we assume Re(z) >> |Im(z)| and Re(z) >> \n 2 — 1 / 4 | (note that Re(z) = dzkjrr, and Im(z) = 
a/A^ — E 2 r/vp 01 ujdt/vf in our integrals), by using the asymptotic form 


we obtain 


Jn (^) 



cos(z — n 7 r /2 — tt/ 4 ) , 


( 73 ) 


x T / n AT cosfz — ( 2 m + 1 — 1 )tt /2 — 7 r/ 4 l 

ij(.z|Re(.z) > 0 ) ~ ttJ-i(z) + 2 zy — £ - 1 2m+ / 7 -^ 


m= — oo 


+oo 


2 ' "7 

= irJ-i(z) + 2 iyJ — sin(z + Ztt/ 2 ^ tt/ 4 ) ^ 

m= — oo 

^ e i(z-n/4+ln/2) = { lP^^z-n/A) 


I l (z\Re(z)<0) = [I- l (-z*)]* ~i l 


Therefore (assuming ujd ^ Too after the integral) 


27 T 


—2T 




where 


oof^) = i b(„,) (£ ; 7 : ) *<„„)*, 


= 


n 2 ir 

Jo 


dpk 

2 i r 


_ ( /-^= c +^+ H- / C-F-) + (s+T + + 5_F_), 


7 A +- B2 7 


A 2 -E 2 


p 2 ir 


A- = 


d<fk . 
——ze 
27r 


C A_ 


( 7 = ■ s + f +- 7 = S_F_) + ( C + F + - C_F_), 


A 2 -E 2 


A 2 -£ 2 


/>27T 


R - / ^ fe R 

B +=L ^ B+ 

~ ~( / C +-P+ H ; : c-F-) - (S+-F+ + s -F-)i 


A 2 . - £ 2 


A 2 - E 2 


B_ = 


p 2 ir 

Jo 


d<Pk ■ i, 


27T 


ie ltfk B_ 


( , = s+F+ - , = s-F_) - (c+F+ - c_F_), 


' ^A 2 . - E 2 yTi - £ 2 

l 


r _ / 
c+ = ' 


( 74 ) 

( 75 ) 

( 76 ) 

( 77 ) 

( 78 ) 

( 79 ) 

( 80 ) 

( 81 ) 

( 82 ) 

( 83 ) 

( 84 ) 

( 85 ) 

( 86 ) 


(87) 
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A+ „ A. 

--C+F+ + 


A \-E 2 


C _ = 


p2n 

JO 


^£l ie i~Pk C _ 

2tt 


A+ „ A_ 

+ -s+F+ - 


A \-E 2 


with 


db) T 

\,p I 


F h = 




k!p r 


: -c F , 

- E 2 

(88) 


(89) 

—s F , 

E 2 

(90) 

e sin[fcpV — 7 t/4], 

(91) 

V F V^V F ' 

(92) 


Now we assume A + = A_ = A, focusing on pure singlet superconductivity, and ignore 0(E 2 / A 2 ) terms (the 
deep-dilute-Shiba-state approximation). Then 


E 

G 0 (E, r) « —t 0 0 a(r) +t x ® a{r) + t z 0 /?(r), 

1 i 

<5(0 = —ji c + F + + c -F- ) - -(s+-F+ - S-F-)(a y cos <p r - a x sin ip r ), 
~ \ % 

P(r) = 2 ^ s+F+ +s ~ F ~) ~ 2^ C+F+ ~ c - F -)( a y cos Fr ~ °x sin <p r ), 


e ~ r /i 


(f = v f /A). 


F \j2 F rv F 

For an FM Shiba lattice, assume the magnetization is along z, we have [cf. Eq. (14)], 

T G o(E, r - r'){-Ja z 0 T 0 )ip(r') = ip(r), 


or, 


E 


— Mi(r - r')ip(r’) = ^ M 0 (r - r')ip(r'), 

r / r' 

Mi(6r = 0) = ^Trjy^J r 0 0<r z , 

Mi(Sr ± 0) = j) r 0 0 [ c + a z - s-(a x cos ip r + a y siny> r )]> 

M 0 (6r = 0) = 1 - ^7rJ T x ®a z , 

M 0 (5r ± 0) = Q j) {-t x 0 [c+<r z - s_ (a x cos <p r + cr,, sin<p r )] 

+ T z 0 [s+<7 z + c_ (a x cos<p r + Gy sin (p r )]}, 


(93) 

(94) 

(95) 

(96) 

(97) 

(98) 

(99) 
( 100 ) 
( 101 ) 

( 102 ) 


p± = k ( p' ) /2nv F , 
c± = c + F + ± c_F_ 

s± = s + F + ± S-F L 




\[2Tvrv F ( 
e ~ r /i 

\/‘2nrv f ( 


(V4 +) cos[4 +) r — 7 t/4] ± \J~k^ F ^ cos[fcp, V — 7 t/4]| , 
( \/ k^ sin[fcp"V — 7 t/4] ± \J~k^ sm[k F V — 7r/4] 1 . 


(103) 

(104) 

(105) 


where 
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In the limit where the Shiba impurities are ultimately dilute, Eq. (98) becomes 


x^( r ) = ( ~ T o ® - T x ® o-Q ) ip{r), 


1 T J O /? 1 

7 r j^p b ~-\^l + -. 


2 vp 

u 

Since we have already assumed the deep Shiba state limit, where 77 ^ 1, we find the two low energy states 


(106) 

(107) 


E± = ± A ( - - 1 ) , ip+ = -j= 


(l\ 

0 

w 


’ *- m T2 


0 \ 
1 
0 

v-v 


Projected into the subspace spanned by these two states, Eq. (98) becomes 

E 


where 


— y]Mi(r - r')^(r') = ^2,M 0 (r - r')i>{r'), 

r' r' 

Mi(5r = 0) = r/a z , 

Mi(6r^0 ) = { l -J)c + a z , 

Mo ( 5 r = 0) = (1 - 77)0-0, 

M 0 (5r ± 0) = (- J) [-c+o-Q + c-(a x cos ip r + a y sin</? r )], 


M i = (V’+,V’-) t M i (^ + ,^_) (i = 0,l), 


(108) 


(109) 

( 110 ) 
( 111 ) 
( 112 ) 

(113) 

(114) 


'ip and a (tildes will be omitted hereafter) are wavefunctions and Pauli matrices under the basis ( < 0 + ,'0_). Fourier 
transforming Eq. (109) we obtain 


Or, 


E 

H e s(q)rp(q) = -£ip(q), 


H eS (q ) = 


e~ iq " r Mi(r) 


-1 


^e-^M 0 (r) 


-Heff(q) = [ 7 ? + d 0 (q)} 1 {[1 - 77 - d 0 (q)\a z + d 1 (q)a x + d 2 (q)a y } , 

d o{q) = (1 A E e_J9 ' r c+(r), 


ry^O 


di(q) = - (^ J )j2 e iq ' r c-{r)s\ 


sin ip r , 


r*/0 


d2{q) = ( ^ J ) E e lqr C-( r ) COS ¥V; 

ry^O 


(115) 

(116) 

(117) 

(118) 

(119) 

( 120 ) 


Note that d\ and are defined with slight differences in the main text in order to make a compact expression. 

We can further construct a Majorana basis by using the fact that and ip- are particle-hole images of each other: 


d M) = ^(i4 + iM = \ 


/1 \ 
1 
1 

V-d 


. 4 =^+-^-) = 2 


/1 \ 
-1 
1 

1 / 


( 121 ) 
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FIG. 5. Chern numbers computed straightforwardly from Hamiltonian (117) as a function of spin-orbit coupling strength, for 
several values of lattice spacing (in units of Fermi wavelength Af) of Shiba impurities. The ill quantization of some of the 
Chern numbers is due to numerical errors. Note that the symbols here are dimensional to be consistent with the main text. 


Under this Major ana basis, the effective Hamiltonian is given by 

H^g\q) = - [n + doiq)} -1 {[1 - q - d 0 (q)}a y + d 2 {q)(J x - di(q)a z } . (122) 

Clearly, when q is one of the inversion symmetric momenta [ISMs, namely, (0, 0), (0,7r), ( 7 r, 0), or ( 7 r, 7 r)], e~ iqr = e iq r , 
therefore d\(q) = ^(g) = 0; the Pfaffian function is given by 


P{q e ISMs) = 1-3—- (123) 

V + d 0 (q ) 

~ do(<jr) (»7 - !)• (124) 

Now return to the effective Hamiltonian Eq. (117). If spin-orbit coupling is absent, then \/q : d\(q) = ^(g) = 0. 
This immediately implies that if there is any sign change in P(q) among the ISMs, the system is gapless with line 
nodes. Rashba spin-orbit coupling potentially lifts these degeneracies. The sign of the Pfaffian formula (124), 
multiplied over both C 4 symmetric momenta is equal to the parity of the Chern number. We have evaluated the 
expression numerically, with the result given in the main text. 

In addition we have evaluated the Chern numbers through numerical integrations by using Hamiltonian (117) 
straightforwardly. Some representative results are shown in Fig. 5. We confirm the existence of high Chern numbers 
associated with this effective Hamiltonian, which complements the phase diagram in terms of the parity of the Chern 
numbers presented in Fig. 3 of the main text. 


Two-dimensional Shiba lattice with helical magnetic order 

In this section we focus on a different type of magnetic order: we assume a two-dimensional helical pattern for the 
magnetic moments. For simplicity we neglect the detailed spatial structure of the Shiba states (which is crucial in the 
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previous section) and investigate the following tight-binding Hamiltonian 


H — T [A n c n (\(Jy)c n + h.c.] + ^ ^ t n s 

n y 5=±l x ,±l y 


c n c n+8 ( ? 


(125) 


where n is a 2D index for the tight-binding sites, c n = (c n ^,c n ^) T is the vector of electron annihilation operators at 
site n, B n stands for the local magnetic moment coupled to the superconductor electron spin, A n stands for the local 
pairing potential, and t n s stands for spin-independent nearest-neighbor hoppings. 

By changing to a rotating basis g n = L^c n with U n defined by 


Ul(B n ■ a)XJ n = B n a z , U}jJ n = cr 0 , Det U n = 1, 

the above Hamiltonian becomes 

H = ^ \ 9n(Bn&z h)9n T n9n [}®y)9n H“ h.C.] Hh ^ ^ tnd 9n^n59n-\-8 


8=±l x ,±l v 


where 


^ n8 — U’^ l U’ n -\-5 5 Det Q n s — 1. 


(126) 


(127) 


(128) 


To further simplify the problem, we set B n = H, A n = A* = A, t n \ x = t* laj = t x , t n \ y = for all n. 

We also assume a certain periodical pattern for the magnetic moments B n . One simplest, nevertheless quite general, 
choice is to let 


n x ,o),8=i x — ^x — 6xp[i(p a ,/2)r a , • <x], 
^‘(n x ,n y ),8=i y =toy = exp[i(p y /2)r y ■ cr], 


(129) 

(130) 


for all n x and n y . Here p x , y stand for rotation angles between neighboring sites and r x , y stand for rotation axes. The 
rest of the ^-matrices are hereby determined by a closed-path formula: 


n { n x ,n y ),8=i x = ^ v) = Vy nv ^ y = exp[i(p x /2)n- n y(r x • <r)fi”»]. 


(131) 


If Cl x and Q y do not commute and p y = 2irp/q [which implies = (— t) p ] with p and q relatively prime integers, 
then Q^ x v ^ = mod q each unit cell consists of q sites. The Bloch Hamiltonian in the Nambu basis is given by 

Hk=g\ 


h(B, k) = 


h(B,k) A \ 

A —h(—B, k)J 9k ' 


(132a) 

y 9oM’ ''' ’ 9q- i 5 fet’ 9q-iM’ ~9o,-kt^ • 

• • 5 9q~ i|— fc-t? 9q— 1 , — &t)) 

(132b) 

h 0 (B,k x ) r(k y ) 0 r{k y ) f \ 

r(fcj / ) t hi(B,k x ) r(k y ) 0 


(132c) 

0 


V T(k y ) 0 r(fc y ) t h q -i(B,k x )J 



= Ba z - p + e ikx + h.c.), 


(132d) 

n y e iky/q . 


(132e) 


The above model is similar, but in general not equivalent, to a model with constant B n and spin-orbit coupling, 
such as the following one: 


hso(B , fe) 


2^2 


h 2 k 

2m 


- p + Bcr z + k x {a. x • a) + kyipLy • cr), 


(133) 


where ol x and oc y are real-valued vectors. To see this, we write the Hamiltonian (133) on a lattice. The hopping 
terms have the form 

{n\h so \n + h= x , y ) = t + b a( • cr. 


(134) 
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FIG. 6. Examples of nodal superconductor band structures with (a) two Fermi points (see Sec. ), (b) four Fermi points (see 
Sec. ). 

If we define 

ti=x,v = + l a */ 2 l 2 > tti=x,y = (t+ lai • cr)/ti, (135) 

where Qi= x , y are unitary and Det Qi = 1, this Hamiltonian formally resembles Hamiltonian (127) except that, impor¬ 
tantly, the fFs such defined do not necessarily satisfy a closed-path equation as Eq. (131). 

Commuting helix rotations 

Obviously, Vt x and Vt y commute when r x = r y = r. In this case 

h(B , k) = Ba z — fi + t x (Q x e lkx + h.c .) + t y (£l y e lky + h.c.). (136) 

Obviously we can rotate the basis again such that r = (sin#, 0, cos#). Then 


7/o u\ (B-v-m -u \ 

k(B,k)-\ _ (B 

(137a) 

m = ii — 2t x cos ^ cos k x — 2t y cos ^ cos k y , 

(137b) 

^ = (2t x sin ^ sin k x + 2t y sin ^ sin k y ) sin #, 

(137c) 

f = (2t x sin ^ sin k x + 2t y sin ^ sin k y ) cos #. 

(137d) 


Let us assume # = 7 t/ 2, the spectrum of the Nambu Hamiltonian is given by 

E = ± \JB 2 + A 2 + rn 2 + u 2 ± 2\JB 2 A 2 + B 2 m 2 + m 2 u 2 . (138) 

The conditions for E = 0 are 

2 t x sin ^ sin k x + 2 t y sin ^ sin k y = 0, (139) 

/i = 2 t x cos ^ cos k x + 2 t y cos ^ cos k y d= \/ B 2 — A 2 . (140) 


Considering B max(A,t) and /i > 0, the spectrum is clearly gapped when |/i — VH 2 — A 2 | > 2£ a ,|cos^ L | + 
2t y | cos otherwise the spectrum is gapless with two Fermi points E(k x o, fc y o) [ see e.g. Fig- 6(a)] that merge when 
the equal sign is taken. 






















































17 


Non-commuting helix rotations 


Now we consider the case when 9l x and 9l y do not commute. The most important feature of the band structure 
is that the subbands created by the multiple sites of a unit cell are not gapped. To see this, let us look at the 
Hamiltonian (132) at k x = 0 or 7r, where e lkx = ±1, and 

g,m) e ifc* + h.c. = ±(gy) + h.c.) = ±2 cos (p x /2). (141) 


This implies that h m (B,k x ) has no dependence on m at k x = 0 or 7r; the unit cell is essentially of a single site. 
Explicitly, the Hamiltonian becomes 


q -1 


Hk x =0,ir(k y ) — ^ dinky 


t (hmi.B, ky ) 


m =0 


Zk /z m ( H, ky) 

q -1 


Qmky 


H / 

~t _ /~t ~f ~ \ ~ 

9 m k = 9m,kV 9m,-kh — 9m,ka = WE e 9 9m',kcr, 

m'= 0 V ^ 

hm(B , fcy) = — p ± 2ic cos(p x /2) + r(27rm + fc^) + r(27rra + fcy)^. 


(142) 

(143) 

(144) 


The crossings of the subbands do not necessarily occur at k y = 0 or 7r unless there is an inversion symmetry. One 
example is that if = exp[i(p 2/ /2)(cr ;E cos cp y + cr y sin^)], then a z Q y a z = fit, and a z h rn (B,k y )cr z = h q - rn (B , ~k y ). 
In other words, if the rotation axis lies in the x-y plane, then the crossings of the subbands occur at k y = 0 or tt. 

As a special case, we set p y = 7r, r y = (1,0,0), = (cos ^, sin^, 0). In this case, each unit cell consists of two 

sites, 9l y = fi^fi^fi^ = fi^, and (we also set t x = t y = t) 


h(B, k) 


f Bcr z - y + t(9l x e lkx + A.c.) 
y — 2 £<j £C sin fcy/2 


—2ta x sin fc y /2 
— p + t(p^e lkx + Ac.) 


(145) 


The spectrum of this Hamiltonian can be solved to be 


E = ±^B 2 + A 2 + m\ + 4 ± ± 2^B 2 A 2 + -B 2 m| + m 2 ^, (146) 

nik = 9 — 2t cos ^ cos k x , (147) 

^4± = 2t(sin ^ sin k x db sin ^). (148) 

The conditions for E = 0 are 

sin ^ sin k x db sin ^ = 0, (149) 

/i = 2 1 cos ^ cos k x db \/ B 2 — A 2 . (150) 


Considering B max(A,t) and p > 0, the spectrum is gapped only when p > |2tcos^| db \/B 2 — A 2 or p < 
— |2tcos db \/B 2 — A 2 . When |/x — VH 2 — A 2 | < |2tcos ^|, there are always four Fermi points [see e.g. Fig. 6(b)]. 

With more general helical patterns, full superconducting gaps can be opened with nontrivial topological properties 
described by Chern numbers [43]. The Chern number, as well as the spectrum for an open-boundary sample, can 
be computed numerically. Unless otherwise stated, we set B = 5, t x = t y = 1, r y = (1,0,0) and r x = (0, —1,0) in 
our computations. We shown in Fig. 7 several examples of the phase diagrams in terms of the Chern numbers with 
different helical patterns. In Fig. 8 we show in addition an example of bulk and edge spectrum for a fixed combination 
of p and A. 


Details of density functional theory analyses 

We performed electronic structure calculations within the density functional theory (DFT) formalism as imple¬ 
mented in the Vienna ab initio simulation package [44]. We used the all-electron projector augmented wave (PAW) 
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Chern number phase diagrams 



FIG. 7. Examples of the Chern number phase diagrams with different helical patterns. 


[45, 46] basis sets with the generalized gradient approximation of Perdew, Burke and Ernzerhof [47] to the exchange 
correlation potential. The Hamiltonian contains scalar relativistic corrections, and the spin-orbit coupling was taken 
into account by the second variation method [48] . 

In this work, we chose the host superconductor to be Pb thin film with a (111) surface, and consider different 
transition metal adatoms. We started by finding the stable configurations of the adatoms on top of the Pb surface. 
To this end we have compared the relaxation energy (per atom) for an extensive collection of possible configurations, 
four of which are shown in Fig. 9. Based on this energetic consideration and for simplicity, we adopted the configuration 
with two adatoms per Pb unit cell (highlighted in Fig. 9) in our following simulations. When the transition metal 
element is chosen to be Fe, the details of the configuration are shown in Fig. 10. In this configuration, two species of 
Fe atoms form a buckled honeycomb structure which gains bonding energy due to the short nearest-neighbor distance. 

To simulate the composite system and consider the effect of Pb, we used 6 layers of Pb atoms as substrate in 
the relaxation calculations with roughly 15 A vacuum space, taking into account spin-orbital coupling and the ferro¬ 
magnetic alignment of the Fe moments. We performed DFT calculations with the stable configuration as shown in 
Fig. 10. Based on the DFT calculations, we constructed the maximally localized Wannier functions [49, 50] for Fe, 
and obtained a tight-binding model with a band structure that agrees well with the DFT result. We then used the 
tight-binding model and added a small 5-wave superconducting pairing term to it. We computed the Chern numbers 
of the thus obtained Bogoliubov-deGennes Hamiltonian to be nonvanishing, as shown in Fig. 4 from the main text. 
For completeness, we present a few more Fermi surfaces with different values of Ep in Fig. 11 to complement Fig. 4 
from the main text. 
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FIG. 8. Low-energy bulk bands (left panel) and dispersion relations for an open-boundary strip (right panel) with p/q = 1/3, 
p x — 27r/3, p — 5.5, A = 0.9. The Chern number is —1 in this example. 



o o 

O Pb 

O o OO 

o o 


© 

o © 
O o 
© 


© 

o © 

© O 
© 


FIG. 9. Relaxation energy of different lattice configurations for different transition metal adatoms. The configuration adopted 
in our simulations is highlighted. 
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FIG. 10. Details of the atomic configuration adopted in our simulations. This configuration has a notable C 3 symmetry. 



lgP(Ef,k) (a.u.) 

FIG. 11. Fermi surfaces with four values of Ep around 0. In the presence of a small (O.OleV) s-wave superconducting pairing 
potential, the Chern number corresponding to the superconducting phase with Ef = 0.18 eV is 3, and the Chern numbers in 
the rest cases are all 1 (see main text Fig. 4). 






























